{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第十四章 复杂问题的随机模拟\n",
    "\n",
    "Code examples from [Think Complexity, 2nd edition](http://greenteapress.com/wp/complexity2), Chapter 7 - 10\n",
    "\n",
    "Copyright 2016 Allen Downey, [MIT License](http://opensource.org/licenses/MIT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这一章在利用上一章提供的基础数据结果和算法，展示如何将这些简单模型应用到复杂的自然和社会问题中去。\n",
    "\n",
    "## 物理模型（physical modeling）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "%matplotlib inline\n",
    "%precision 3\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import thinkplot\n",
    "\n",
    "from matplotlib import rc\n",
    "rc('animation', html='html5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 扩散\n",
    "\n",
    "很多自然现象都可以用反应-扩散（Reaction-Diffusion）模型来描述。物质间相互作用并改变彼此的浓度是反应，而物质的浓度在空间上的分布随时间的变化则可以用扩散来描述。对于扩散这种最基本的现象，我们用在最简单的格子空间中，可以用计算一个单元和周围单元的差分（difference）来实现。每个时间步，实际工作就是计算差分，乘以一个常数，加回矩阵。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Cell2D import Cell2D, Cell2DViewer\n",
    "from scipy.signal import correlate2d\n",
    "\n",
    "class Diffusion(Cell2D):   # 继承自Cell2D\n",
    "    \"\"\"Diffusion Cellular Automaton.\"\"\"\n",
    "\n",
    "    kernel = np.array([[0, 1, 0],\n",
    "                       [1,-4, 1],\n",
    "                       [0, 1, 0]])\n",
    "    \n",
    "    def __init__(self, n, m=None, r=0.1):\n",
    "        \"\"\"Initializes the attributes.\n",
    "\n",
    "        n: number of rows\n",
    "        m: number of columns\n",
    "        r: diffusion rate constant\n",
    "        \"\"\"\n",
    "        self.r = r\n",
    "        m = n if m is None else m\n",
    "        self.array = np.zeros((n, m), np.float)\n",
    "\n",
    "    def step(self):\n",
    "        \"\"\"Executes one time step.\"\"\"\n",
    "        c = correlate2d(self.array, self.kernel, mode='same')\n",
    "        self.array += self.r * c   # 乘比例 r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "来看看这个kernel的实际含义。假设中心单元坐标为$x_i, y_j$，单元的值记为$u_{i,j}:=u(x_i, y_j)$，则其上下左右的邻居为$x_i, y_{j+1}$，$x_i, y_{j-1}$，$x_{i-1},y_j$和$x_{i+1},y_j$；其相应的单元值为$u_{i,j+1}$，$u_{i,j-1}$，$u_{i-1}$，$u_{i+1}$。那么对应kernel\n",
    "$$\n",
    "\\left[\n",
    "\\begin{array}{ccc}\n",
    "0&1&0\\\\\n",
    "1&-4&1\\\\\n",
    "0&1&0\n",
    "\\end{array}\n",
    "\\right],\n",
    "$$\n",
    "correlate2d的实际作用是\n",
    "$$\n",
    "u_{i,j}=u_{i,j}+(-4u_{i,j}+u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1})*r.\n",
    "$$\n",
    "而这个从另一个方向来看，我们知道热传导方程：\n",
    "$$\n",
    "\\frac{\\partial u}{\\partial t} = -c \\nabla^2 u = -c \\left(\\frac{\\partial^2 u}{\\partial x^2} + \\frac{\\partial^2u}{\\partial y^2}\\right),\n",
    "$$\n",
    "其中$c>0$为扩散率。对其做五点差分离散，并记网格宽度为$h$，时间步长为$\\Delta t:=t_{n+1}-t_n$，$u^{(n)}_{i, j}:=u(x_i, y_j, t_n)$，则其差分离散格式为：\n",
    "$$\n",
    "\\frac{u^{(n+1)}_{i,j}-u^{(n)}_{i,j}}{\\Delta t}=c\\left(\\frac{u^{(n)}_{i+1,j}+u^{(n)}_{i-1,j}-2u^{(n)}_{i,j}}{h^2}+\\frac{u^{(n)}_{i,j+1}+u^{(n)}_{i,j-1}-2u^{(n)}_{i,j}}{h^2}\\right)\n",
    "$$\n",
    "即\n",
    "$$\n",
    "u^{(n+1)}_{i,j}=u^{(n)}_{i,j}+\\frac{\\Delta t c}{h^2}\\left(-4u^{(n)}_{i,j}+u^{(n)}_{i+1,j}+u^{(n)}_{i-1,j}+u^{(n)}_{i,j+1}+u^{(n)}_{i,j-1}\\right),\n",
    "$$\n",
    "所以实际上是对$r=\\frac{\\Delta t c}{h^2}$的五点差分离散求解。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class DiffusionViewer(Cell2DViewer):\n",
    "    cmap = plt.get_cmap('Reds')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们的计算区域的中心先放置一块物质，然后来模拟它的扩散情况。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff = Diffusion(n=9, r=0.1)   # 9行9列，扩撒率r包含了时间和空间离散尺度\n",
    "diff.add_cells(3, 3, '111', '111', '111')   # 初始状态\n",
    "viewer = DiffusionViewer(diff)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "初值状态："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anim = viewer.animate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们可以看到物质是如何扩散的。注意这里单元中的值不再保持0或1，而是连续的浮点数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 分别画出初值，第5步和第10步。\n",
    "diff = Diffusion(n=9, r=0.1)\n",
    "diff.add_cells(3, 3, '111', '111', '111')\n",
    "viewer = DiffusionViewer(diff)\n",
    "\n",
    "thinkplot.preplot(cols=3)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "viewer.step(5)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(3)\n",
    "viewer.step(10)\n",
    "viewer.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "大家注意，一直扩散下去，比如扩散无穷步，解会是什么？这要看边界如何设置。也就是物质扩散到边界了，如何和区域之外进行交换。在这里，collerate2D实质上是把边界之外看作全0的。因此，这里等价地施加了全0边界条件。也就是说，若$t\\to \\infty$，整个区域都会趋于0。一个模拟在$t\\to \\infty$的解（若存在），称为稳态。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff = Diffusion(n=9, r=0.1)\n",
    "diff.add_cells(3, 3, '111', '111', '111')\n",
    "viewer = DiffusionViewer(diff)\n",
    "\n",
    "thinkplot.preplot(cols=3)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "viewer.step(5)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(3)\n",
    "viewer.step(100)   # 测试一下稳态\n",
    "viewer.draw()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 反应扩散（reaction-diffusion）\n",
    "\n",
    "现在我们考虑两种物质的扩散以及之间的相互作用。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def island(a, val, noise=None):\n",
    "    \"\"\"Adds an island in the middle of the array.\n",
    "    \n",
    "    增加一块物质。\n",
    "            \n",
    "    val: height of the island\n",
    "    \n",
    "    val表示物质的浓度。现在单元中不再是0|1值，可以是浮点值。\n",
    "    \n",
    "    noise: magnitude of random noise\n",
    "    \n",
    "    noise是背景噪音强度，或者扰动。模拟物质的自然浓度起伏。\n",
    "    \"\"\"\n",
    "    noise = val if noise is None else noise\n",
    "    n, m = a.shape\n",
    "    r = min(n, m) // 20\n",
    "    a[n//2-r:n//2+r, m//2-r:m//2+r] = val\n",
    "    a += noise * np.random.random((n, m)) \n",
    "    # 加上一个噪音起伏. 这个为了增强非线性效应，更接近物理实际。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里我们考虑反应扩散模型。是两种物质间的交互。用两个数值矩阵表示。\n",
    "\n",
    "参考 [Sims](http://www.karlsims.com/rd.html)。\n",
    "\n",
    "每一步的反应规则如下。物质 $A$ 和物质 $B$ 分别以 $r_a$ 和 $r_b$ 的扩散率扩散，然后 $A$ 的浓度每个时间步以比率 $f$ 关于 $1-A$ 增长($A$ 越小，增长越快，$A$ 越接近 $1$，增长越慢，因此 $A$ 是一个 $0$，$1$之间的数)，$B$ 的浓度以比率 $k+f$ 关于 $B$ 下降，注意 $B$ 浓度下降率总是高于 $A$ 的上升率的。然后 $A$ 和 $B$ 之间还有一个复杂的非线性效应 $AB^2$，对于 $A$ 是下降，对于 $B$ 是上升。它们反应了 $A$ 和 $B$ 两种物质之间的化学反应变化。$A$ 和 $B$ 在全空间的浓度分布变化整体上就受到这一组合的扩散和反应两种效应的影响。\n",
    "\n",
    "$\\Delta A = r_a \\nabla^2 A - AB^2 + f (1-A) $\n",
    "\n",
    "$\\Delta B = r_b \\nabla^2 B + AB^2 - (k+f) B $\n",
    "\n",
    "这里 $\\nabla^2$ 是 Laplace 算子。这个模型的具体导出，参见图灵的论文。不同的扩散率、上升和下降率会对模型的结果造成影响，就是下面的参数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ReactionDiffusion(Cell2D):\n",
    "    \"\"\"Reaction-Diffusion Cellular Automaton.\"\"\"\n",
    "\n",
    "    kernel = np.array([[.05, .2, .05],\n",
    "                       [ .2, -1, .2],\n",
    "                       [.05, .2, .05]])\n",
    "    \n",
    "    options = dict(mode='same', boundary='wrap')\n",
    "\n",
    "    def __init__(self, n, m=None, params=(0.5, 0.25, 0.035, 0.057)):\n",
    "        \"\"\"Initializes the attributes.\n",
    "\n",
    "        n: number of rows\n",
    "        m: number of columns\n",
    "        params: tuple of (Da, Db, f, k)\n",
    "        \"\"\"\n",
    "        \n",
    "        self.params = params\n",
    "        m = n if m is None else m\n",
    "        self.array = np.ones((n, m), dtype=float)\n",
    "\n",
    "        self.array2 = np.zeros((n, m), dtype=float)\n",
    "        island(self.array2, val=0.1, noise=0.1)\n",
    "        \n",
    "    def step(self):\n",
    "        \"\"\"Executes one time step.\"\"\"\n",
    "        A = self.array\n",
    "        B = self.array2\n",
    "        ra, rb, f, k = self.params\n",
    "        \n",
    "        cA = correlate2d(A, self.kernel, **self.options)\n",
    "        cB = correlate2d(B, self.kernel, **self.options)\n",
    "        reaction = A * B**2\n",
    "        self.array += ra * cA - reaction + f * (1-A) \n",
    "        self.array2 += rb * cB + reaction - (f+k) * B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由于不再是 0|1 值，所以绘图的时候不能用两色了，而要搞个 0 ~ 1 插值反应浓度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class RDViewer(Cell2DViewer):\n",
    "    \"\"\"Generates images and animations.\"\"\"\n",
    "    \n",
    "    cmapu = plt.get_cmap('Reds')\n",
    "    cmapv = plt.get_cmap('Blues')\n",
    "\n",
    "    options = dict(alpha=0.7,\n",
    "                  interpolation='bicubic')\n",
    "        \n",
    "    def __init__(self, viewee):\n",
    "        \"\"\"Initializes the attributes.\n",
    "        \n",
    "        viewee: the object to be represented\n",
    "        \"\"\"\n",
    "        self.viewee = viewee\n",
    "        self.imu = None\n",
    "        self.imv = None\n",
    "        self.hlines = None\n",
    "        self.vlines = None\n",
    "\n",
    "    def draw(self, grid=False):\n",
    "        \"\"\"Draws the cells.\"\"\"\n",
    "        au = self.viewee.array.copy()\n",
    "        av = self.viewee.array2.copy()\n",
    "        \n",
    "        n, m = av.shape\n",
    "        plt.axis([0, m, 0, n])\n",
    "        plt.xticks([])\n",
    "        plt.yticks([])\n",
    "\n",
    "        self.options['extent'] = [0, m, 0, n]\n",
    "        self.imu = plt.imshow(au, cmap=self.cmapu, **self.options)\n",
    "        self.imv = plt.imshow(av, cmap=self.cmapv, **self.options)\n",
    "\n",
    "    def animate_func(self, i):\n",
    "        \"\"\"Draws one frame of the animation.\"\"\"\n",
    "        if i > 0:\n",
    "            self.step(iters=100)\n",
    "\n",
    "        self.imu.set_array(self.viewee.array)\n",
    "        self.imv.set_array(self.viewee.array2)\n",
    "        return (self.imu, self.imv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面给了三组参数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params1 = 0.5, 0.25, 0.035, 0.057   # white spots\n",
    "params2 = 0.5, 0.25, 0.055, 0.062   # coral\n",
    "params3 = 0.5, 0.25, 0.039, 0.065   # blue spots\n",
    "\n",
    "rd = ReactionDiffusion(100, params=params3)\n",
    "viewer = RDViewer(rd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这是一个以B为主的初值，A是随机的噪音，B初始有很多杂点，以及中心一个岛。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anim = viewer.animate(frames=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "会如何演化？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面秀一下其他参数的图像。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_viewer(f, k, n=100):\n",
    "    \"\"\"Makes a ReactionDiffusion viewer with given parameters.\n",
    "    \"\"\"\n",
    "    params = 0.5, 0.25, f, k\n",
    "    rd = ReactionDiffusion(n, params=params)\n",
    "    viewer = RDViewer(rd)\n",
    "    return viewer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "蓝色背景上绘制粉色的条纹和点。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viewer = make_viewer(0.035, 0.057)\n",
    "\n",
    "thinkplot.preplot(cols=3)\n",
    "viewer.step(1000)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "viewer.step(2000)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(3)\n",
    "viewer.step(4000)\n",
    "viewer.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "或粉色背景上的蓝条纹。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viewer = make_viewer(0.055, 0.062)\n",
    "\n",
    "thinkplot.preplot(cols=3)\n",
    "viewer.step(1000)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "viewer.step(2000)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(3)\n",
    "viewer.step(4000)\n",
    "viewer.draw()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "粉色背景的蓝色点"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viewer = make_viewer(0.039, 0.065)\n",
    "\n",
    "thinkplot.preplot(cols=3)\n",
    "viewer.step(1000)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "viewer.step(2000)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(3)\n",
    "viewer.step(4000)\n",
    "viewer.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 渗透（Percolation）\n",
    "\n",
    "现在我们来模拟多孔介质，比如防水瓦，的渗透。每个单元有概率 $p$ 为多孔性的。每个时间步，一个单元会不会变湿，需要满足：1. 自身是多孔的；2. 至少有一个邻居是湿的。我们研究对一个具体的 $p$，水能渗透到底层的概率是多少。\n",
    "\n",
    "多孔的单元用状态1表达，湿单元用状态5。所以若至少有一个邻居单元是湿的，那么邻居状态之和应大于5。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import correlate2d\n",
    "\n",
    "class Percolation(Cell2D):\n",
    "    \"\"\"Percolation Cellular Automaton.\"\"\"\n",
    "\n",
    "    kernel = np.array([[0, 1, 0],\n",
    "                       [1, 0, 1],\n",
    "                       [0, 1, 0]])\n",
    "    # 这里可以看出邻居指的还是上下左右单元。\n",
    "    \n",
    "    def __init__(self, n, m=None, p=0.5, seed=None):\n",
    "        \"\"\"Initializes the attributes.\n",
    "\n",
    "        n: number of rows\n",
    "        m: number of columns\n",
    "        p: probability of porousness\n",
    "        \"\"\"\n",
    "        self.p = p\n",
    "        m = n if m is None else m\n",
    "        if seed is not None:\n",
    "            np.random.seed(seed)\n",
    "        self.array = np.random.choice([0, 1], (n, m), p=[1-p, p])\n",
    "        # 依概率来均匀分布多孔介质。 \n",
    "        \n",
    "        # fill the top row with wet cells\n",
    "        # 最顶上的一排单元，设置成湿的。（这是一个边界条件）\n",
    "        self.array[0] = 5\n",
    "\n",
    "    def step(self):\n",
    "        \"\"\"Executes one time step.\"\"\"\n",
    "        a = self.array\n",
    "        c = correlate2d(a, self.kernel, mode='same')\n",
    "        # 自己是多孔的（a==1），邻居是湿的（c>=5），则变湿。\n",
    "        self.array[(a==1) & (c>=5)] = 5\n",
    "        \n",
    "    def num_wet(self, cols=None):\n",
    "        \"\"\"Total number of wet cells.\n",
    "\n",
    "        统计湿单元总数。\n",
    "        \n",
    "        cols: number of columns to select\n",
    "        \"\"\"\n",
    "        a = self.array[:, :cols]\n",
    "        return np.sum(a == 5)\n",
    "    \n",
    "    def num_porous(self, cols=None):\n",
    "        \"\"\"Total number of porous cells.\n",
    "            \n",
    "        统计多孔单元总数。\n",
    "        \"\"\"\n",
    "        a = self.array[:, :cols]\n",
    "        return np.sum(a == 1)\n",
    "    \n",
    "    def bottom_row_wet(self):\n",
    "        \"\"\"Number of wet cells in the bottom row.\n",
    "        \n",
    "        统计最后一排单元中的湿单元数。\n",
    "        \"\"\"\n",
    "        return np.sum(self.array[-1] == 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "显示部分和之前没什么区别："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class PercolationViewer(Cell2DViewer):\n",
    "    \"\"\"Draws and animates a Percolation object.\"\"\"\n",
    "    cmap = plt.get_cmap('Blues')\n",
    "    options = dict(alpha=0.6,\n",
    "                   interpolation='none', \n",
    "                   vmin=0, vmax=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "看3步。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = 22\n",
    "perc = Percolation(10, p=0.5, seed=seed)\n",
    "viewer = PercolationViewer(perc)\n",
    "\n",
    "thinkplot.preplot(cols=3)\n",
    "viewer.step()\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "viewer.step()\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(3)\n",
    "viewer.step()\n",
    "viewer.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面的函数用这个模拟来做模拟，如果水最终能到达底部，返回True，否则False。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_perc(perc):\n",
    "    \"\"\"Run a percolation model.\n",
    "    \n",
    "    Runs until water gets to the bottom row or nothing changes.\n",
    "    \n",
    "    returns: (boolean, number of steps)\n",
    "    \"\"\"\n",
    "    num_wet = perc.num_wet()\n",
    "\n",
    "    num_steps = 0\n",
    "    while True:\n",
    "        perc.step()\n",
    "        num_steps += 1\n",
    "\n",
    "        # 如果最后一行湿了，退出。\n",
    "        if perc.bottom_row_wet():\n",
    "            return True, num_steps\n",
    "        # 如果到了某一步湿格子的总数不再增加，表明渗透停止。\n",
    "        new_num_wet = perc.num_wet()\n",
    "        if new_num_wet == num_wet:\n",
    "            return False, num_steps\n",
    "\n",
    "        num_wet = new_num_wet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "0.5的渗透率，11步渗透完毕。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "perc = Percolation(10, p=0.5, seed=seed)\n",
    "flag, num_steps = test_perc(perc)\n",
    "print(flag, num_steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "动态过程："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "perc = Percolation(10, p=0.5, seed=seed)\n",
    "viewer = PercolationViewer(perc)\n",
    "anim = viewer.animate(frames=num_steps+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对给定的 $p$，我们可以计算最终透水的概率："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimate_prob_percolating(p=0.5, n=100, iters=100):\n",
    "    \"\"\"Estimates the probability of percolating.\n",
    "    \n",
    "    估计完全渗透的概率。\n",
    "    \n",
    "    p: probability that a cell is permeable\n",
    "    n: int number of rows and columns\n",
    "    iters: number of arrays to test\n",
    "    \n",
    "    returns: float probability\n",
    "    \"\"\"\n",
    "    count = 0\n",
    "    for i in range(iters):\n",
    "        perc = Percolation(n, p=p)\n",
    "        flag, _ = test_perc(perc)\n",
    "        if flag:\n",
    "            count += 1\n",
    "        \n",
    "    return count / iters   # 总共测试了多少次，渗透了多少次。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对 $p=0.55$ ，最终透水是多少？（默认铺了 100 层）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fraction = estimate_prob_percolating(p=0.55)\n",
    "print(fraction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在假设你是一个建筑师，打算用防水瓦来构建防水层。如果高质量的瓦很昂贵（$p$ 值很低），你能否用质量不太高的瓦（$p$ 值很大），通过多铺几层来达到防水的效果？比如刚才我们发现 $p = 0.5$（一半的瓦是坏的），实际的透水率并不高。我们进一步试一下 $p = 0.6$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fraction = estimate_prob_percolating(p=0.6)\n",
    "print(fraction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "呵呵，和我们想象的有点不一样，再试试 $p = 0.65$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fraction = estimate_prob_percolating(p=0.65)\n",
    "print(fraction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以上结果说明，在 0.6 附近有一个阈值，会使得透水率从接近零极速上升到接近1，我们用一个随机游走实验把它找出来："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_critical(p=0.6, n=100, iters=100):\n",
    "    \"\"\"Estimate p_crit by random walk.\n",
    "    \n",
    "    returns: list of p that should wander around p_crit\n",
    "    \"\"\"\n",
    "    ps = [p]\n",
    "    for i in range(iters):\n",
    "        perc = Percolation(n=n, p=p)\n",
    "        flag, _ = test_perc(perc)   # 这是第二项输出不要了的意思么\n",
    "        if flag:   # 如果有渗透，则减小概率\n",
    "            p -= 0.005\n",
    "        else:  # 否则增加概率  \n",
    "            p += 0.005\n",
    "        ps.append(p)\n",
    "    return ps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们先取 $n = 50$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time ps = find_critical(n=50, iters=1000)\n",
    "thinkplot.plot(ps)\n",
    "np.mean(ps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "也就是说在 0.59 附近，有一个快速转变的阈值。我们试一下 $n = 100$，以验证是否可以通过多铺几层将强防水性。这里我们将尝试次数降低，以免等得太久。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time ps = find_critical(n=100, iters=200)\n",
    "thinkplot.plot(ps)\n",
    "np.mean(ps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "提高 $n$ 似乎并不影响阀值。我们再来试一个 $n = 200$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time ps = find_critical(n=200, iters=40)\n",
    "thinkplot.plot(ps)\n",
    "np.mean(ps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们现在可以得到结论：最佳的质量应该控制在比 $p = 0.59$ 更优，而如果单片瓦片质量不足 $p = 0.59$，提高铺设层数并不能提高最终工程质量。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 智能体建模（Agent Based Models）\n",
    "\n",
    "之前的模型，都可以称为基于规则的(rule-based)建模，主要特点是基于简单的规则。这一章的智能体建模(agent based models)的主要应用是对人类行为进行建模。模型的对象会从外部收取信息，作出判定，然后给出反应。这些智能体可以是空间上的分布或一个网络，并且局部相邻的个体间会有交互。但交互的信息往往是残缺的和局部的。不像之前的模型，个体都是相同的，这里的个体间往往有差异。并且智能模型在整体上和个体上，都往往会包含随机因素。\n",
    "从1970年代以来，智能体建模变成经济和社会科学领域的一种重要建模工具。特别是对一些没有非平衡状态的动态系统(当然它对平衡系统也有应用)，以及帮助我们了解个体决定和整体行为之间的关系。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "%matplotlib inline\n",
    "%precision 3\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import animation\n",
    "\n",
    "import thinkplot\n",
    "from thinkstats2 import Cdf\n",
    "from thinkstats2 import RandomSeed\n",
    "\n",
    "from matplotlib import rc\n",
    "rc('animation', html='html5')\n",
    "\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Schelling's model\n",
    "\n",
    "早期的重要工作之一是1971年Thomas Schelling发表的Dynamic Models of Segregation。这个工作提供了一个种族隔离现象的简单模型。他的模型世界是一个正方形网格，每个网格代表一个房子。然后其间生活的人群分两种，分别用红和蓝表示，大致上具有相同的数量，每个人占据一个房子，使网格便红或蓝。然后大约有10%的房子是空置状态。\n",
    "\n",
    "在任何时刻，每个人都有两个状态，happy或unhappy，由他们的邻居决定。每个格子空间上有8个邻居。在一种版本的模型中，每个个体如果有指定个数以上同色的邻居，则会处于happy状态；若是不到指定个数，则会unhappy。\n",
    "\n",
    "整个计算区域先随机初始化，然后根据每个个体的happy与否的状态来决定变化，happy的个体，不做任何变化，unhappy的个体，则随机地移动到一个空置网格。显然，这种模型会导致一些分离现象，但这种分离的程度仍然出乎一般的意料。随着模拟，很快就出现一些相似个体构成的小团体，然后这些小团体会增长和合并，知道最终形成少量的几个大型的团体，于是大部分个体，最终都和相同的邻居居住。\n",
    "\n",
    "如果你不了解这个过程，并只看到了结果，你可能会下结论这是一个种族主义盛行的社会。但实际上的原因可能只是每个人都想和自己类似的人做邻居(是个体选择，而不是一种宏观控制造成的。)因为每个人都不想在邻居间显得突兀。当然，实际的人类行为原比此模型复杂，因此反而难以给出这种简单的描述。\n",
    "\n",
    "种族主义是一种复杂的人类问题，很难想像这么简单的模型可以阐明这个问题。但它给出了一个基于整体和个体之间关系的重要论据：如果你在一个城市观察到种族隔离现象，你不能判定是否是个体的种族主义造成的，以及这个城市的居民是否是种族主义者。\n",
    "\n",
    "我们这里也要强调Schelling模型的局限：该模型给出了种族隔离的一种可能的原因，但并没有分析种族隔离的实际原因。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Cell2D import Cell2D, Cell2DViewer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where 函数的功能是返回满足条件的坐标："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def locs_where(condition):\n",
    "    \"\"\"Find cells where a logical array is True.\n",
    "    \n",
    "    condition: logical array\n",
    "    \n",
    "    returns: array with one set of coordinates per row\n",
    "    \"\"\"\n",
    "    # transpose是转置\n",
    "    return np.transpose(np.nonzero(condition))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 这段代码给出了locs_where的测试结果\n",
    "# 返回全部>0.5的A的坐标\n",
    "A=(np.random.random(9)).reshape(3,3)\n",
    "R=locs_where(A>0.5)\n",
    "print(A)\n",
    "print(R)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们还是基于Cell2D来构建Schelling类。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import correlate2d\n",
    "\n",
    "class Schelling(Cell2D):\n",
    "    \"\"\"Represents a grid of Schelling agents.\"\"\"\n",
    "    # 这里我们用周期边界条件\n",
    "    options = dict(mode='same', boundary='wrap')\n",
    "\n",
    "    kernel = np.array([[1, 1, 1],\n",
    "                       [1, 0, 1],\n",
    "                       [1, 1, 1]], dtype=np.int8)\n",
    "    \n",
    "    def __init__(self, n, m=None, p=0.5):\n",
    "        \"\"\"Initializes the attributes.\n",
    "\n",
    "        n: number of rows\n",
    "        m: number of columns\n",
    "        p: threshold on the fraction of similar neighbors\n",
    "        \n",
    "        这里p是happy与否的阀值。0.5表示如果有一半以上的邻居是和自己类似的，\n",
    "        则happy。\n",
    "        \"\"\"\n",
    "        self.p = p\n",
    "        m = n if m is None else m\n",
    "        # 0 is empty, 1 is red, 2 is blue\n",
    "        choices = [0, 1, 2]\n",
    "        # 各自的比例\n",
    "        probs = [0.1, 0.45, 0.45]\n",
    "        # 数据初始化\n",
    "        self.array = np.random.choice(choices, (n, m), p=probs).astype(np.int8)\n",
    "\n",
    "    def count_neighbors(self):\n",
    "        \"\"\"Surveys neighboring cells.\n",
    "        \n",
    "        returns: tuple of (occupied, frac_red, frac_same)\n",
    "        \n",
    "        occupied: logical array indicating occupied cells\n",
    "        frac_red: fraction of red neighbors around each cell\n",
    "        frac_same: for each cell, the fraction of similar neighbors\n",
    "        \"\"\"\n",
    "        a = self.array\n",
    "        \n",
    "        red = a==1   # 逻辑判断，返回一个和a一样维数的矩阵, ==1成立的地方是1，\n",
    "        blue = a==2   # 其余为0。其他两行类似。\n",
    "        occupied = a!=0\n",
    "\n",
    "        # 分别计算每个位置的红色邻居和全部邻居\n",
    "        num_red = correlate2d(red, self.kernel, **self.options)\n",
    "        num_neighbors = correlate2d(occupied, self.kernel, **self.options)\n",
    "\n",
    "        # 计算和邻居的相似率\n",
    "        frac_red = num_red / num_neighbors\n",
    "        frac_blue = 1 - frac_red\n",
    "        # 红色的地方填充红色的邻居比，蓝色的地方填充蓝色。\n",
    "        frac_same = np.where(red, frac_red, frac_blue)\n",
    "        \n",
    "        # 没邻居的设成unhappy\n",
    "        frac_same[num_neighbors == 0] = 0\n",
    "        frac_red[num_neighbors == 0] = 0\n",
    "        \n",
    "        return occupied, frac_red, frac_same\n",
    "\n",
    "    def segregation(self):\n",
    "        \"\"\"Computes the average fraction of similar neighbors.\n",
    "        \n",
    "        计算平均邻居相似度，这个可以作为种族隔离程度的一个指标。\n",
    "        \n",
    "        returns: fraction of similar neighbors, averaged over cells\n",
    "        \"\"\"\n",
    "        occupied, _, frac_same = self.count_neighbors()\n",
    "        return np.sum(frac_same) / np.sum(occupied)\n",
    "    \n",
    "    def step(self):\n",
    "        \"\"\"Executes one time step.\n",
    "        \n",
    "        模型发展一步。\n",
    "                \n",
    "        returns: fraction of similar neighbors, averaged over cells\n",
    "        \"\"\"\n",
    "        a = self.array\n",
    "        \n",
    "        # 定位unhappy的格子\n",
    "        occupied, _, frac_same = self.count_neighbors()\n",
    "        # 有人住，且相同邻居比小于p，则unhappy\n",
    "        unhappy_locs = locs_where(occupied & (frac_same < self.p))\n",
    "\n",
    "        # 定位空格子\n",
    "        empty = a==0\n",
    "        num_empty = np.sum(empty)\n",
    "        empty_locs = locs_where(empty)\n",
    "\n",
    "        # 乱序unhappy的格子，使整个迁移过程更均匀\n",
    "        if len(unhappy_locs):\n",
    "            np.random.shuffle(unhappy_locs)\n",
    "            \n",
    "        # 随机选一个目的地\n",
    "        for source in unhappy_locs:\n",
    "            # 随便找一个，不用相邻\n",
    "            i = np.random.randint(len(empty_locs))\n",
    "            dest = tuple(empty_locs[i])\n",
    "\n",
    "            # 移动\n",
    "            a[dest] = a[tuple(source)]\n",
    "            a[tuple(source)] = 0\n",
    "            empty_locs[i] = source\n",
    "        \n",
    "        num_empty2 = np.sum(a==0)\n",
    "        assert num_empty == num_empty2\n",
    "        \n",
    "        # 返回隔离度\n",
    "        return np.sum(frac_same) / np.sum(occupied)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`make_cmap` 构建颜色对映矩阵。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.colors import LinearSegmentedColormap\n",
    "\n",
    "def make_cmap(color_dict, vmax=None, name='mycmap'):\n",
    "    \"\"\"Makes a custom color map.\n",
    "    \n",
    "    color map就是给每一个数值对应一个颜色。可以离散，可以连续。\n",
    "    这个函数的功能就是将一组已经规定了(数值,颜色值)对应关系的color_dict\n",
    "    转换成python能接受的color map。\n",
    "    color_dict: map from numbers to colors\n",
    "    vmax: high end of the range, 最大的接受数值\n",
    "    name: string name for map\n",
    "    \n",
    "    If vmax is None, uses the max value from color_dict\n",
    "    \n",
    "    returns: pyplot color map\n",
    "    \"\"\"\n",
    "    if vmax is None:\n",
    "        vmax = max(color_dict.keys())\n",
    "    # 关于vmax做个归一化，避免数值过于奇异\n",
    "    colors = [(value/vmax, color) for value, color in color_dict.items()]\n",
    "    \n",
    "    # 实际上做个分段线性插值，并转换为matplotlib能接受的cmap格式\n",
    "    cmap = LinearSegmentedColormap.from_list(name, colors)\n",
    "    \n",
    "    return cmap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "基本上除了一个cmap，SchellingViewer就是Cell2DViewer。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SchellingViewer(Cell2DViewer):\n",
    "    # colors from http://colorbrewer2.org/#type=qualitative&scheme=Accent&n=5\n",
    "    colors = ['#7fc97f','#beaed4','#fdc086','#ffff99','#386cb0']\n",
    "    cmap = make_cmap({0:'white', 1:colors[2], 2:colors[4]})\n",
    "    options = dict(interpolation='none', alpha=0.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "看个例子："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = Schelling(n=10)\n",
    "viewer = SchellingViewer(grid)\n",
    "viewer.draw()\n",
    "grid.segregation()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对一个更大规模的例子，给出一个动态变化："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = Schelling(n=100)\n",
    "viewer = SchellingViewer(grid)\n",
    "anim = viewer.animate(frames=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "合并的程度和速度，都出乎意料。我们观察一下 2 到 10 步之间的隔离度变化。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RandomSeed(17)\n",
    "\n",
    "grid = Schelling(n=100, p=0.3)\n",
    "print(grid.segregation())\n",
    "\n",
    "thinkplot.preplot(cols=3)\n",
    "viewer = SchellingViewer(grid)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "for i in range(2):\n",
    "    viewer.step()\n",
    "print(grid.segregation())\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(3)\n",
    "for i in range(8):\n",
    "    viewer.step()\n",
    "print(grid.segregation())\n",
    "viewer.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对不同的p做测试。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RandomSeed(17)\n",
    "\n",
    "for p in [0.5, 0.4, 0.3, 0.2]:\n",
    "    grid = Schelling(n=100, p=p)\n",
    "    segs = [grid.step() for i in range(20)]\n",
    "    thinkplot.plot(segs, label='p = %.1f' % p)\n",
    "    print(p, segs[-1], segs[-1] - p)\n",
    "    \n",
    "thinkplot.config(xlabel='Time steps', ylabel='Segregation',\n",
    "                loc='lower right', ylim=[0, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以看到，对$p = 0.3$，实际上邻居相似度的接受程度是很宽松的，有30%的类似邻居就行了，即便如此，最后稳定下来的整体隔离度仍然高达80%。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 糖域（Sugarscape）\n",
    "\n",
    "糖域是一种针对社会行为的人工建模。本节的内容主要来自Joshua Epstein和Robert Axtell在1996年写的书Growing Artifical Societies.\n",
    "\n",
    "糖域是一种基于2维网格的智能体建模，智能体的目标是采集和积累sugar，这是一种对社会财富的抽象。一些网格会产生更多的sugar，而一些智能体具有更好的发现财富的能力。这种模型主要用于分析社会财富的分布，特别是为何会有不平等的倾向。\n",
    "\n",
    "每一个网格都有一个容量(capacity)，代表了能容纳sugar的最大量。在初始状态，有两个高sugar区，都拥有4单位sugar，然后依次被3,2,1包围。然后我们随机抛下400个智能体，每个智能体有三个被随机分配的属性：\n",
    "\n",
    "1. Sugar: 每个智能体被初始随机分配了5～25单位的sugar；\n",
    "2. 消耗(Metabolism)：每个时间步，每个智能体会消耗的sugar，在1~4之间；\n",
    "3. 视野(Vision)：每个智能体能看到周围一定范围内的单元的sugar分布，然后移动到最多的单元。这个范围对不同的智能体是1～6随机分配。\n",
    "\n",
    "每个时间步，每个智能体以随机的乱序进行观察和移动。并遵守下面规则：\n",
    "\n",
    "1. 智能体在4个坐标方向上各观察$k$个单元，$k$是视野范围；\n",
    "2. 选择还未被占据的拥有最多sugar的单元，如果这样的单元不止一个，则选择最近的一个。如果距离也相等，则随机选其中一个。\n",
    "3. 移动到那个单元，并采集那里的sugar，积累到自己的sugar总量里，并留下一个空单元。\n",
    "4. 每个智能体根据自己的消耗值减少sugar，如果总量小于0,则被饿死，并删除。\n",
    "\n",
    "每执行完若干步，每个单元会增加一些sugar，通常是1单位，但不超过每个单元的容量。\n",
    "\n",
    "每个智能体都朝着sugar最多的区域移动。有大视距的智能体移动最快，而小视距的智能体则会在次级平台游走，直到偶尔离更高层近到能看见。\n",
    "\n",
    "那些出生在低sugar区域的智能体更容易饿死，除非他们有很高的视距或者初始sugar。而拥有高消耗或短视距的智能体，更容易饿死。\n",
    "\n",
    "当每个时间步，若每个单元增加一个单位的sugar，这不足以支撑400个初始智能体的消耗。于是智能体总量会快速下降，然后会缓和一些，最终稳定在250左右。\n",
    "\n",
    "100个时间步后，更多的好运智能体，也就是出生在高sugar地区并拥有地消耗的智能体，更容易活下来。而活到这个时候的智能体，几乎可以长生不死。并积累无穷的财富。\n",
    "\n",
    "这个模型是一个简单的经济学模型，并提供了一些参数之间的关系。比如智能体的变化率、属性、稳定的智能体数量之间的关系。结论是自然选择的一个例子，适应性(fitness)更好的个体，更容易存活。\n",
    "\n",
    "模型也揭示了一些财富不平衡现象。某些智能体总是比其他能更快收集财富。但很难将财富的分布有什么特异，因为它没有稳态，总是在不停地变化。\n",
    "\n",
    "然而，如果我们给智能体有限的寿命，则模型会提供稳定的财富分布。然后我们可以观察诸参数对这种分布的影响。\n",
    "\n",
    "在这种版本的模型下，每个智能体都会有一个随时间步增长的年龄，但年龄到达寿命上限，是随机分配的60-100的范围，智能体会老死。不管智能体是饿死还是老死的，它会被一个重新随机分配属性的新个体替代，因此智能体总数是守恒的。\n",
    "\n",
    "从250个初始财富差不多的个体开始，运行500步，并分别在100、200、300、400步时绘制sugar的分布，则收敛的趋势是明显的。25%的个体的财富值是10，中位是差不多20，而但位于75%高位的个体的财富差不多是40，而最高个体拥有150。在log分布下，财富的分布更像一个右偏的高斯或正态分布。事实上，在真实的世界，财富的分布是log下的正态分布，也就是长尾分布的。\n",
    "\n",
    "尽管糖域用于解释财富不平衡的原因还是不够的，但从中观察到的财富不平衡至少说明了大部分经济体，哪怕非常简单，也会不平等，而且这种不平等是极难消除或减轻的。\n",
    "\n",
    "显然糖域是复杂的多的模型。相比之前而言。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_locs(n, m):\n",
    "    \"\"\"Makes array where each row is an index in an `n` by `m` grid.\n",
    "    \n",
    "    产生n行m列的矩阵下标顺序列表\n",
    "    \n",
    "    n: int number of rows\n",
    "    m: int number of cols\n",
    "    \n",
    "    returns: NumPy array\n",
    "    \"\"\"\n",
    "    left = np.repeat(np.arange(m), n)\n",
    "    right = np.tile(np.arange(n), m)\n",
    "    return np.transpose([left, right])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "make_locs(3, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_visible_locs(vision):\n",
    "    \"\"\"Computes the kernel of visible cells.\n",
    "    \n",
    "    计算可视范围。这里故意使观察方向随机混乱。\n",
    "    \n",
    "    vision: int distance\n",
    "    \"\"\"\n",
    "    def make_array(d):\n",
    "        \"\"\"Generates visible cells with increasing distance.\"\"\"\n",
    "        a = np.array([[-d, 0], [d, 0], [0, -d], [0, d]])\n",
    "        np.random.shuffle(a)\n",
    "        return a\n",
    "                     \n",
    "    arrays = [make_array(d) for d in range(1, vision+1)]\n",
    "    return np.vstack(arrays)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "糖域的实现:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Sugarscape(Cell2D):\n",
    "    \"\"\"\n",
    "    Represents an Epstein-Axtell Sugarscape.\n",
    "    \n",
    "    表现的是 EA 的糖域模型。\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, n, **params):\n",
    "        \"\"\"Initializes the attributes.\n",
    "\n",
    "        属性初始化，n是行列数，其余参数全部塞进字典params。\n",
    "        n: number of rows and columns\n",
    "        params: dictionary of parameters\n",
    "        \"\"\"\n",
    "        self.n = n\n",
    "        self.params = params\n",
    "        \n",
    "        # track variables\n",
    "        # 统计智能体总数\n",
    "        self.agent_count_seq = []\n",
    "    \n",
    "        # make the capacity array\n",
    "        # 规定单元糖容量\n",
    "        self.capacity = self.make_capacity()\n",
    "        \n",
    "        # initially all cells are at capacity\n",
    "        # 初始糖量就是初始容量\n",
    "        self.array = self.capacity.copy()\n",
    "        \n",
    "        # make the agents\n",
    "        # 产生智能体\n",
    "        self.make_agents()\n",
    "        \n",
    "    def make_capacity(self):\n",
    "        \"\"\"\n",
    "        Makes the capacity array.\n",
    "        \n",
    "        产生容量矩阵。\n",
    "        \"\"\"\n",
    "        \n",
    "        # compute the distance of each cell from the peaks. \n",
    "        # 人为构建了两个峰，(15,15)和(35,35)\n",
    "        # x = [1,2,...,n-1]\n",
    "        x = np.arange(self.n)\n",
    "        # 这个和Matlab中的meshgrid差不多，产生二维的整点坐标。\n",
    "        rows, cols = np.meshgrid(x, x, indexing='ij')\n",
    "        \n",
    "        # each cell in `rows` contains its own `i` coordinate\n",
    "        # each cell in `cols` contains its `j` coordinate\n",
    "        # 分别放在rows和cols中，这个就是算欧氏距离的吧。\n",
    "        dist1 = np.hypot(rows-15, cols-15)\n",
    "        dist2 = np.hypot(rows-35, cols-35)\n",
    "        \n",
    "        # each cell in `dist` contains its distance to the closer peak\n",
    "        # 取更近的一个作为距离\n",
    "        dist = np.minimum(dist1, dist2)\n",
    "        \n",
    "        # cells in the capacity array are set according to dist from peak\n",
    "        # a 先全部取0，然后距离<21，取1，<16，取2，<11，取3，<6，取4。\n",
    "        a = np.zeros((self.n, self.n), np.float)\n",
    "        a[dist<21] = 1\n",
    "        a[dist<16] = 2\n",
    "        a[dist<11] = 3\n",
    "        a[dist<6] = 4\n",
    "\n",
    "        return a\n",
    "        \n",
    "    def make_agents(self):\n",
    "        \"\"\"Makes the agents.\"\"\"\n",
    "        \n",
    "        # determine where the agents start and generate locations\n",
    "        # n,m取参数中的starting_box，或者取自己矩阵array的shape，也就是\n",
    "        # 分布在全糖域\n",
    "        n, m = self.params.get('starting_box', self.array.shape)\n",
    "        # 域中每个点的坐标\n",
    "        locs = make_locs(n, m)\n",
    "        # 随机混乱一下\n",
    "        np.random.shuffle(locs)\n",
    "\n",
    "        # make the agents\n",
    "        # 智能体初始总数取参数中的num_agents或400\n",
    "        num_agents = self.params.get('num_agents', 400)\n",
    "        # 其初始坐标从locs中一一随机抽取而成\n",
    "        self.agents = [Agent(locs[i], self.params) \n",
    "                       for i in range(num_agents)]\n",
    "        \n",
    "        # keep track of which cells are occupied\n",
    "        # 标记一下哪个坐标已经被占领\n",
    "        self.occupied = set(agent.loc for agent in self.agents)\n",
    "            \n",
    "    def grow(self):\n",
    "        \"\"\"\n",
    "        Adds sugar to all cells and caps them by capacity.\n",
    "        全部单元格增长糖分，受限于每个格子的上限。\n",
    "        \"\"\"\n",
    "        # 增长率（每个时间步每单元各增长的糖量）设为参数中的grow_rate或1\n",
    "        grow_rate = self.params.get('grow_rate', 1)\n",
    "        # 每个格子增长糖分，然后受限于各自容量\n",
    "        self.array = np.minimum(self.array + grow_rate, self.capacity)\n",
    "        \n",
    "    def look_around(self, center, vision):\n",
    "        \"\"\"\n",
    "        寻找附近最多糖分的单元（这个操作明显应该放在Agent里）\n",
    "        Finds the visible cell with the most sugar.\n",
    "        中心出发的单元\n",
    "        center: tuple, coordinates of the center cell\n",
    "        视距\n",
    "        vision: int, maximum visible distance\n",
    "        返回附近糖分最多单元的坐标\n",
    "        returns: tuple, coordinates of best cell\n",
    "        \"\"\"\n",
    "        # find all visible cells\n",
    "        # 以(0,0)为中心，构建可视范围\n",
    "        locs = make_visible_locs(vision)\n",
    "        # 偏移到实际的中心，再和边界做一个周期处理\n",
    "        locs = (locs + center) % self.n\n",
    "        \n",
    "        # convert rows of the array to tuples\n",
    "        # 将全体可视格子转成坐标\n",
    "        locs = [tuple(loc) for loc in locs]\n",
    "        \n",
    "        # select unoccupied cells\n",
    "        # 挑出全部未被占领的格子\n",
    "        empty_locs = [loc for loc in locs if loc not in self.occupied]\n",
    "        \n",
    "        # if all visible cells are occupied, stay put\n",
    "        # 如果视距内全部格子都被占领，停着不动。\n",
    "        if len(empty_locs) == 0:\n",
    "            return center\n",
    "        \n",
    "        # look up the sugar level in each cell\n",
    "        # 看一下空余格子的糖分\n",
    "        \n",
    "        t = [self.array[loc] for loc in empty_locs]\n",
    "        \n",
    "        # find the best one and return it\n",
    "        # (in case of tie, argmax returns the first, which\n",
    "        # is the closest)\n",
    "        # 挑出含糖最多的格子。如果有多个格子含糖相同，则挑出最近的一个。\n",
    "        # 这里最近是通过argmax的特性，总是返回序列中标号小的一个。同时，\n",
    "        # make_visible_locs里面是按距离从近到远排序的。（同样距离是乱序）\n",
    "        i = np.argmax(t)\n",
    "        return empty_locs[i]\n",
    "    \n",
    "    def harvest(self, loc):\n",
    "        \"\"\"Removes and returns the sugar from `loc`.\n",
    "        删除被占据格子的糖分。\n",
    "        \n",
    "        loc: tuple coordinates\n",
    "        \"\"\"\n",
    "        sugar = self.array[loc]\n",
    "        self.array[loc] = 0\n",
    "        return sugar\n",
    "    \n",
    "    def step(self):\n",
    "        \"\"\"\n",
    "        Executes one time step.\n",
    "        \"\"\"\n",
    "        # repalce属性代表了智能体死亡了是不是要增加一个\n",
    "        # 默认是False，如果是True，则智能体总数是守恒的。\n",
    "        replace = self.params.get('replace', False)\n",
    "        \n",
    "        # loop through the agents in random order\n",
    "        # 全部智能体随机排个序\n",
    "        random_order = np.random.permutation(self.agents)\n",
    "        # 对每一个智能体\n",
    "        for agent in random_order:\n",
    "            # mark the current cell unoccupied\n",
    "            # 把当前的单元标为未占领\n",
    "            self.occupied.remove(agent.loc)\n",
    "            \n",
    "            # execute one step\n",
    "            # 每个智能体也发展一步\n",
    "            agent.step(self)\n",
    "\n",
    "            # if the agent is dead, remove from the list\n",
    "            # 如果智能体死亡了，则删除（饿死或老死）\n",
    "            if agent.is_starving() or agent.is_old():\n",
    "                self.agents.remove(agent)\n",
    "                # 如果replace属性是True，死一个就要加一个\n",
    "                if replace:\n",
    "                    self.add_agent()\n",
    "            else:\n",
    "                # otherwise mark its cell occupied\n",
    "                # 否则就把新的位置占据了。\n",
    "                self.occupied.add(agent.loc)\n",
    "\n",
    "        # update the time series\n",
    "        # 更新智能体数量变化的随时间变化的序列\n",
    "        self.agent_count_seq.append(len(self.agents))\n",
    "        \n",
    "        # grow back some sugar\n",
    "        # 普遍增加糖\n",
    "        self.grow()\n",
    "        # 返回当前智能体总数\n",
    "        return len(self.agents)\n",
    "    \n",
    "    def add_agent(self):\n",
    "        \"\"\"Generates a new random agent.\n",
    "        增加一个随机的智能体。\n",
    "        returns: new Agent\n",
    "        \"\"\"\n",
    "        # 产生一个随机位置和随机属性的智能体\n",
    "        new_agent = Agent(self.random_loc(), self.params)\n",
    "        # 增加进智能体序列\n",
    "        self.agents.append(new_agent)\n",
    "        # 标记为已占据\n",
    "        self.occupied.add(new_agent.loc)\n",
    "        return new_agent\n",
    "    \n",
    "    def random_loc(self):\n",
    "        \"\"\"Choose a random unoccupied cell.\n",
    "        随机挑选一个未被占据的单元。\n",
    "        returns: tuple coordinates\n",
    "        \"\"\"\n",
    "        while True:\n",
    "            # 在糖域范围不断随机抽取\n",
    "            loc = tuple(np.random.randint(self.n, size=2))\n",
    "            # 如果未必占据，则返回坐标\n",
    "            if loc not in self.occupied:\n",
    "                return loc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SugarscapeViewer(Cell2DViewer):\n",
    "    \"\"\"Generates visualization and animation of Sugarscape.\"\"\"\n",
    "    \n",
    "    cmap = plt.get_cmap('YlOrRd')\n",
    "    \n",
    "    options = dict(interpolation='none', alpha=0.8,\n",
    "                   vmin=0, vmax=9)\n",
    "    \n",
    "    def draw(self, grid=False):\n",
    "        \"\"\"Draws the array and any other elements.\n",
    "        \n",
    "        grid: boolean, whether to draw grid lines\n",
    "        \"\"\"\n",
    "        self.draw_array(self.viewee.array, origin='lower')\n",
    "        self.draw_agents()\n",
    "\n",
    "    def draw_agents(self):\n",
    "        \"\"\"Plots the agents.\n",
    "        \"\"\"\n",
    "        xs, ys = self.get_coords()\n",
    "        self.points = plt.plot(xs, ys, '.', color='red')[0]\n",
    "        \n",
    "    def animate_func(self, i):\n",
    "        \"\"\"Draws one frame of the animation.\"\"\"\n",
    "        Cell2DViewer.animate_func(self, i)\n",
    "        xs, ys = self.get_coords()\n",
    "        self.points.set_data(np.array([xs, ys]))\n",
    "        return self.im, self.points\n",
    "    \n",
    "    def get_coords(self):\n",
    "        \"\"\"Gets the coordinates of the agents.\n",
    "        \n",
    "        Transforms from (row, col) to (x, y).\n",
    "        \n",
    "        returns: tuple of sequences, (xs, ys)\n",
    "        \"\"\"\n",
    "        agents = self.viewee.agents\n",
    "        rows, cols = np.transpose([agent.loc for agent in agents])\n",
    "        xs = cols + 0.5\n",
    "        ys = rows + 0.5\n",
    "        return xs, ys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "智能体也是一个类。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Agent:\n",
    "    \n",
    "    def __init__(self, loc, params):\n",
    "        \"\"\"Creates a new agent at the given location.\n",
    "        在指定位置产生一个随机属性的智能体\n",
    "        \n",
    "        loc: tuple coordinates\n",
    "        params: dictionary of parameters\n",
    "        \"\"\"\n",
    "        self.loc = tuple(loc)\n",
    "        # 年龄设为0\n",
    "        self.age = 0\n",
    "\n",
    "        # extract the parameters\n",
    "        # 导入诸变量或采用默认值\n",
    "        max_vision = params.get('max_vision', 6)\n",
    "        max_metabolism = params.get('max_metabolism', 4)\n",
    "        min_lifespan = params.get('min_lifespan', 10000)\n",
    "        max_lifespan = params.get('max_lifespan', 10000)\n",
    "        min_sugar = params.get('min_sugar', 5)\n",
    "        max_sugar = params.get('max_sugar', 25)\n",
    "        \n",
    "        # choose attributes\n",
    "        # 实际设置具体智能体的属性参数，都是随机生成\n",
    "        self.vision = np.random.random_integers(max_vision)\n",
    "        # 注意消耗、寿命和初始糖是浮点数\n",
    "        self.metabolism = np.random.uniform(1, max_metabolism)\n",
    "        self.lifespan = np.random.uniform(min_lifespan, max_lifespan)\n",
    "        self.sugar = np.random.uniform(min_sugar, max_sugar)\n",
    "\n",
    "    def step(self, env):\n",
    "        \"\"\"Look around, move, and harvest.\n",
    "        这里外部环境（糖域）的参数通过env读入\n",
    "        env: Sugarscape\n",
    "        \"\"\"\n",
    "        # 这几句话显得略蠢，一开始的逻辑没组织好就会这样。\n",
    "        self.loc = env.look_around(self.loc, self.vision)\n",
    "        # 收入是新格子采集来的糖，支出是自身的消耗\n",
    "        self.sugar += env.harvest(self.loc) - self.metabolism\n",
    "        # 年龄增加\n",
    "        self.age += 1\n",
    "\n",
    "    def is_starving(self):\n",
    "        \"\"\"\n",
    "        检查饿死了么。\n",
    "        Checks if sugar has gone negative.\n",
    "        \"\"\"\n",
    "        return self.sugar < 0\n",
    "    \n",
    "    def is_old(self):\n",
    "        \"\"\"\n",
    "        Checks if lifespan is exceeded.\n",
    "        检查老死了么。\n",
    "        \"\"\"\n",
    "        return self.age > self.lifespan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "做个测试，50x50的格子上，放入400个智能体。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# env是一个糖域实体\n",
    "env = Sugarscape(50, num_agents=400)\n",
    "# 看一下初始状态\n",
    "viewer = SugarscapeViewer(env)\n",
    "viewer.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "全体智能体的视距均匀分布在1到6之间。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cdf = Cdf(agent.vision for agent in env.agents)\n",
    "thinkplot.Cdf(cdf)\n",
    "\n",
    "# 有图为证\n",
    "cdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "也可以看一下消耗的分布。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cdf = Cdf(agent.metabolism for agent in env.agents)\n",
    "thinkplot.Cdf(cdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "初始携糖量也是5-25之间的均匀分布。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cdf = Cdf(agent.sugar for agent in env.agents)\n",
    "thinkplot.Cdf(cdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "然后整个糖域的动态发展如下。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anim = viewer.animate(frames=500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "整个糖域能容纳的智能体总数是一个定量。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thinkplot.plot(env.agent_count_seq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "看看2到10步的动态变化。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RandomSeed(17)\n",
    "\n",
    "env = Sugarscape(50, num_agents=400)\n",
    "viewer = SugarscapeViewer(env)\n",
    "\n",
    "thinkplot.preplot(cols=3)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "for i in range(2):\n",
    "    viewer.step()\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(3)\n",
    "for i in range(98):\n",
    "    viewer.step()\n",
    "viewer.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "然而，不同的出身和天赋，会导致财富不均么？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cdf = Cdf(agent.sugar for agent in env.agents)\n",
    "thinkplot.Cdf(cdf)\n",
    "thinkplot.config(xlabel='Wealth', ylabel='CDF')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 有限寿命的糖域 （Sugarscape with finite lifespans）\n",
    "\n",
    "这一次设置糖域有250个智能体，寿命在60-100之间，每个死去的智能体总会被新生代替，因此智能体总数不变。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "env = Sugarscape(50, num_agents=250,\n",
    "                 min_lifespan=60, max_lifespan=100, \n",
    "                 replace=True)\n",
    "\n",
    "viewer = SugarscapeViewer(env)\n",
    "anim = viewer.animate(frames=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "100步以后，财富的分布开始不均匀，大部分智能体只有少量的财富，而少数则有很多。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cdf = Cdf(agent.sugar for agent in env.agents)\n",
    "thinkplot.Cdf(cdf)\n",
    "thinkplot.config(xlabel='Wealth', ylabel='CDF')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "再来一次，但每100步记录一下财富分布，直到500步。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RandomSeed(17)\n",
    "\n",
    "env = Sugarscape(50, num_agents=250,\n",
    "                 min_lifespan=60, max_lifespan=100, \n",
    "                 replace=True)\n",
    "\n",
    "cdfs = []\n",
    "for i in range(5):\n",
    "    [env.step() for i in range(100)]\n",
    "    cdf = Cdf(agent.sugar for agent in env.agents)\n",
    "    cdfs.append(cdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "大概200步起，财富分布开始稳定。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thinkplot.preplot(cols=2)\n",
    "\n",
    "thinkplot.Cdfs(cdfs[:-1], color='gray', alpha=0.3)\n",
    "thinkplot.Cdf(cdfs[-1])\n",
    "thinkplot.config(xlabel='Weath', ylabel='CDF')\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "thinkplot.Cdfs(cdfs[:-1], color='gray', alpha=0.3)\n",
    "thinkplot.Cdf(cdfs[-1])\n",
    "thinkplot.config(xlabel='Weath', ylabel='CDF', xscale='log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
